1. Introdução

O objetivo deste relatório é documentar todas as análises e etapas de modelagem que levaram ao desenvolvimento do protótipo do AgriPredict+, com dados públicos fornecidos pela Fundação MS (FMS).

2. Análise Exploratória e Tratamento de Dados

2.1. Precipitação Pluviométrica

Os dados brutos fornecidos pela FMS já haviam sofrido um tratamento preliminar por parte do Rodrigo Stelzer, o cientista de dados com quem se iniciou o projeto. Dele, tive acesso a um arquivo csv com dados de precipitação pluviométrica, a saber: clima_dia.csv

A primeira etapa da análise desse arquivo consiste em carregá-lo e verificar sua estrutura:

clima_dia <- read.csv("data-raw/clima_dia.csv",
                      stringsAsFactors = F) # arquivo do Stelzer
str(clima_dia)
'data.frame':   4018 obs. of  10 variables:
 $ date                 : chr  "2008-01-01" "2008-01-02" "2008-01-03" "2008-01-04" ...
 $ AmambaÃ.             : num  NA NA NA NA NA NA NA NA NA NA ...
 $ Anaurilandia         : num  0 0 39 0 0 0 0 0 0 0 ...
 $ Campo_Grande         : num  0 0 3.4 2.2 4.4 28.6 0.4 29.6 0 1.2 ...
 $ Dourados             : num  0 0.2 1.6 28.6 12.8 0 0 0 0 0 ...
 $ Ivinhema             : num  8.8 44.4 0 5.4 0 0 0 0 0 0 ...
 $ Maracaju             : num  0.4 0 0 30.6 1.6 0.4 0.2 0 0 0 ...
 $ Rio_Brilhante        : num  NA NA NA NA NA NA NA NA NA NA ...
 $ SÃ.o_Gabriel_do_Oeste: num  19 7.6 0.4 5 0 1 0.2 30 0.8 9.4 ...
 $ SidrolÃ.ndia         : num  NA NA NA NA NA NA NA NA NA NA ...
summary(clima_dia)
     date              AmambaÃ.        Anaurilandia     Campo_Grande    
 Length:4018        Min.   :  0.000   Min.   :  0.00   Min.   :  0.000  
 Class :character   1st Qu.:  0.000   1st Qu.: 15.00   1st Qu.:  0.000  
 Mode  :character   Median :  0.000   Median : 22.50   Median :  0.000  
                    Mean   :  4.001   Mean   : 30.51   Mean   :  4.053  
                    3rd Qu.:  1.600   3rd Qu.: 40.00   3rd Qu.:  1.800  
                    Max.   :147.600   Max.   :180.00   Max.   :111.000  
                    NA's   :449       NA's   :3556     NA's   :308      
    Dourados          Ivinhema          Maracaju      Rio_Brilhante   
 Min.   :  0.000   Min.   :  0.000   Min.   : 0.000   Min.   :  0.00  
 1st Qu.:  0.000   1st Qu.:  0.000   1st Qu.: 0.000   1st Qu.:  0.00  
 Median :  0.000   Median :  0.000   Median : 0.000   Median :  0.00  
 Mean   :  3.471   Mean   :  3.567   Mean   : 2.828   Mean   :  3.99  
 3rd Qu.:  1.000   3rd Qu.:  0.600   3rd Qu.: 0.600   3rd Qu.:  1.20  
 Max.   :186.000   Max.   :146.000   Max.   :83.400   Max.   :135.00  
 NA's   :407       NA's   :372       NA's   :286      NA's   :369     
 SÃ.o_Gabriel_do_Oeste  SidrolÃ.ndia    
 Min.   : 0.000        Min.   :  0.000  
 1st Qu.: 0.000        1st Qu.:  0.000  
 Median : 0.000        Median :  0.000  
 Mean   : 1.988        Mean   :  3.247  
 3rd Qu.: 0.400        3rd Qu.:  0.800  
 Max.   :95.000        Max.   :128.800  
 NA's   :398           NA's   :927

Nota-se o elevado número de NAs para todos os locais, notadamente Anaurilândia, e que o campo de data está representado como character, de modo que alguns ajustes se fazem então necessários:

clima_dia$date <- as.Date(clima_dia$date)
round(colSums(is.na(clima_dia))/nrow(clima_dia) * 100, 0) # % NAs
                 date              AmambaÃ.          Anaurilandia 
                    0                    11                    89 
         Campo_Grande              Dourados              Ivinhema 
                    8                    10                     9 
             Maracaju         Rio_Brilhante SÃ.o_Gabriel_do_Oeste 
                    7                     9                    10 
         SidrolÃ.ndia 
                   23 
fNames <- c("Date", "Amambai", "Anaurilandia", "Campo.Grande", "Dourados",
            "Ivinhema", "Maracaju", "Rio.Brilhante", "Sao.Grabiel.do.Oeste",
            "Sidrolandia")
names(clima_dia) <- fNames
clima_dia$Anaurilandia <- NULL # me livro dessa cidade, pois tem muitos NAs
clima_dia <- na.omit(clima_dia)
summary(clima_dia)
      Date               Amambai         Campo.Grande    
 Min.   :2008-10-01   Min.   :  0.000   Min.   :  0.000  
 1st Qu.:2010-10-30   1st Qu.:  0.000   1st Qu.:  0.000  
 Median :2013-06-11   Median :  0.000   Median :  0.000  
 Mean   :2013-06-25   Mean   :  4.111   Mean   :  4.174  
 3rd Qu.:2015-10-02   3rd Qu.:  1.400   3rd Qu.:  2.400  
 Max.   :2018-07-25   Max.   :147.600   Max.   :111.000  
    Dourados          Ivinhema          Maracaju      Rio.Brilhante    
 Min.   :  0.000   Min.   :  0.000   Min.   : 0.000   Min.   :  0.000  
 1st Qu.:  0.000   1st Qu.:  0.000   1st Qu.: 0.000   1st Qu.:  0.000  
 Median :  0.000   Median :  0.000   Median : 0.000   Median :  0.000  
 Mean   :  3.141   Mean   :  3.818   Mean   : 2.728   Mean   :  3.992  
 3rd Qu.:  0.950   3rd Qu.:  0.800   3rd Qu.: 0.400   3rd Qu.:  1.200  
 Max.   :186.000   Max.   :146.000   Max.   :81.400   Max.   :135.000  
 Sao.Grabiel.do.Oeste  Sidrolandia     
 Min.   : 0.000       Min.   :  0.000  
 1st Qu.: 0.000       1st Qu.:  0.000  
 Median : 0.000       Median :  0.000  
 Mean   : 2.436       Mean   :  3.363  
 3rd Qu.: 0.600       3rd Qu.:  0.600  
 Max.   :95.000       Max.   :128.800

A fim de agregar os dados trimestralmentre, teremos de fazer uso do pacote xts:

clima.dia.xts <- xts(clima_dia[, -1], order.by = clima_dia$Date)
clima.qrt.xts <- apply.quarterly(clima.dia.xts, mean)
ts_plot(clima.qrt.xts, type = "single")
periodicity(clima.qrt.xts)
Quarterly periodicity from 2008-12-31 to 2018-07-25 
frequency(clima.qrt.xts)
[1] 1

Do gráfico anterior e de entrevistas com o Sr. Alex Melotto, presidente da FMS, parece factível agregar as precipitações dos diversos locais em uma única variável. Um mapa de correlação nos ajudará a avaliar melhor essa hipótese:

clima.cor <- cor(clima.qrt.xts)
corrplot.mixed(clima.cor)

ts_plot(clima.qrt.xts[, c("Rio.Brilhante", "Ivinhema")])

À exceção de São Gabriel do Oeste, fica evidente a lucidez dessa proposta, i.e. agregar as precipitações dos locais, que é o que se faz a seguir:

clima.qrt <- as.data.frame(clima.qrt.xts)
clima.qrt$Date <- as.Date(rownames(clima.qrt))
rownames(clima.qrt) <- NULL
clima.qrt$quarter <- quarter(clima.qrt$Date)
clima.qrt$harvest <- year(clima.qrt$Date)

clima.qrt <- clima.qrt %>% select(harvest, quarter, everything(), -Date)
clima.qrt$rain.mm <- rowMeans(clima.qrt[, 3:ncol(clima.qrt)])
clima.qrt <- clima.qrt %>% select(harvest, quarter, rain.mm)
str(clima.qrt)
'data.frame':   36 obs. of  3 variables:
 $ harvest: num  2008 2009 2009 2009 2009 ...
 $ quarter: int  4 1 2 3 4 1 2 3 4 4 ...
 $ rain.mm: num  3.55 5.23 1.49 2.94 5.2 ...

Finalmente, por solicitação do próprio Sr. Alex Melotto, seria melhor que os dados trimestrais de precipitação fossem traduzidos em classes, e.g. “sem chuva”, “pouca chuva”, “chuva média”, “muita chuva”, em vez de números.

# Níveis de precipitação
rain.level <- cut2(clima.qrt$rain.mm, g = 4)
levels(rain.level) <- c("norain", "low", "med", "high")
clima.qrt$rain.level <- as.character(rain.level)
clima.qrt$rain <- paste0('q', clima.qrt$quarter, '.', clima.qrt$rain.level)
clima.qrt$rain <- factor(clima.qrt$rain, levels = c(
  'q1.norain', 'q1.low', 'q1.med', 'q1.high',
  'q2.norain', 'q2.low', 'q2.med', 'q2.high',
  'q3.norain', 'q3.low', 'q3.med', 'q3.high',
  'q4.norain', 'q4.low', 'q4.med', 'q4.high'))

# Um pouco de feature engineering...
clima.qrt$rain.level <- NULL
dummies <- dummy_cols(clima.qrt)
dummies <- dummies %>% select(harvest, 
                              rain_q1.low, rain_q1.med, rain_q1.high,
                              rain_q2.norain, rain_q2.low, rain_q2.med,
                              rain_q3.norain, rain_q3.low,
                              rain_q4.low, rain_q4.med, rain_q4.high)

clima.qrt.final <- dummies %>% 
  group_by(harvest) %>% summarise_if(is.numeric, max)
str(clima.qrt.final)
Classes 'tbl_df', 'tbl' and 'data.frame':   11 obs. of  12 variables:
 $ harvest       : num  2008 2009 2010 2011 2012 ...
 $ rain_q1.low   : int  0 0 0 0 0 0 0 1 0 0 ...
 $ rain_q1.med   : int  0 0 0 0 1 0 0 0 1 1 ...
 $ rain_q1.high  : int  0 1 1 0 0 1 1 0 0 0 ...
 $ rain_q2.norain: int  0 1 1 0 0 0 0 0 1 0 ...
 $ rain_q2.low   : int  0 0 0 0 0 1 1 0 0 1 ...
 $ rain_q2.med   : int  0 0 0 0 1 0 0 0 0 0 ...
 $ rain_q3.norain: int  0 0 1 0 1 1 0 0 0 1 ...
 $ rain_q3.low   : int  0 1 0 0 0 0 1 1 1 0 ...
 $ rain_q4.low   : int  0 0 1 0 0 0 0 0 0 0 ...
 $ rain_q4.med   : int  1 0 0 0 0 1 1 1 1 0 ...
 $ rain_q4.high  : int  0 1 0 1 1 0 0 0 0 1 ...

Finalmente, salvo o arquivo resultante:

write_csv(clima.qrt.final, "clima.csv")

2.2. Produtividade das culturas de soja e milho

Os dados brutos de produtividade das culturas de soja já haviam sofrido um tratamento preliminar por parte do Rodrigo Stelzer, o cientista de dados com quem se iniciou o projeto. Dele, tive acesso a um arquivo csv com esses dados, a saber: crop.csv

A primeira etapa da análise desse arquivo consiste em carregá-lo e verificar sua estrutura:

crop <- read.csv("data-raw/crop.csv", encoding="UTF-8", 
                 stringsAsFactors = F) # arquivo do Stelzer
str(crop)
'data.frame':   14436 obs. of  6 variables:
 $ id          : int  0 1 2 3 4 5 6 7 8 9 ...
 $ material    : chr  "fts campo mourão rr" "na 5909 rg" "brs 282" "brs 284" ...
 $ harvest     : int  2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
 $ site        : chr  "antônio joão" "antônio joão" "antônio joão" "antônio joão" ...
 $ productivity: num  66 65 63.6 63.5 59.9 58.4 57.8 57.3 56.6 56.2 ...
 $ species     : chr  "soy" "soy" "soy" "soy" ...
summary(crop)
       id         material            harvest         site          
 Min.   :   0   Length:14436       Min.   :2008   Length:14436      
 1st Qu.:1804   Class :character   1st Qu.:2010   Class :character  
 Median :3608   Mode  :character   Median :2013   Mode  :character  
 Mean   :3612                      Mean   :2013                     
 3rd Qu.:5413                      3rd Qu.:2015                     
 Max.   :7444                      Max.   :2018                     
  productivity      species         
 Min.   :  2.60   Length:14436      
 1st Qu.: 58.60   Class :character  
 Median : 72.80   Mode  :character  
 Mean   : 84.41                     
 3rd Qu.:114.10                     
 Max.   :198.60                     
length(unique(crop$material))
[1] 835
barplot(table(crop$material))

Da análise anterior, nota-se um grande número de sementes distintas, cujas classes estão nitidamente desbalanceadas, o que pode ser um problema para o treinamento dos algoritmos de machine learning.

A seguir, mais algumas análises e gráficos de exploração:

crop$harvest <- factor(crop$harvest)
crop$id <- NULL

boxplot(productivity ~ harvest, data = subset(crop, species == 'soy'),
        main = 'Produtividade da Soja')

boxplot(productivity ~ harvest, data = subset(crop, species == 'corn'),
        main = 'Produtividade do Milho')

table(crop$harvest)/nrow(crop)*100

     2008      2009      2010      2011      2012      2013      2014 
 1.780272 11.277362 12.835966 10.092823  7.460515 10.944860 11.360488 
     2015      2016      2017      2018 
 9.843447 12.524245  9.434746  2.445276 
cumsum(table(crop$harvest)/nrow(crop)*100) # até 2015 para treinamento, resto para teste
      2008       2009       2010       2011       2012       2013 
  1.780272  13.057634  25.893599  35.986423  43.446938  54.391798 
      2014       2015       2016       2017       2018 
 65.752286  75.595733  88.119978  97.554724 100.000000 
g <- ggplot(data = subset(crop, species == 'soy'), 
            aes(x = harvest, y = productivity)) + geom_boxplot() + 
  facet_wrap(~ site, ncol = 4) + labs(title = 'Produtividade da Soja')
g + theme(axis.text.x = element_text(angle = 90, hjust = 1))

g1 <- ggplot(data = subset(crop, species == 'corn'), 
            aes(x = harvest, y = productivity)) + geom_boxplot() + 
  facet_wrap(~ site, ncol = 4) + labs(title = 'Produtividade do Milho')
g1 + theme(axis.text.x = element_text(angle = 90, hjust = 1))

crop.sum <- crop %>% group_by(species, site) %>% 
  summarise(avgProd = mean(productivity), sdProd = sd(productivity)) %>%
  arrange(desc(avgProd))
crop.sum
# A tibble: 32 x 4
# Groups:   species [2]
   species site                 avgProd sdProd
   <chr>   <chr>                  <dbl>  <dbl>
 1 corn    sidrolândia            128.    26.6
 2 corn    bonito                 126.    16.7
 3 corn    amambai                118.    26.2
 4 corn    são gabriel do oeste   116.    22.7
 5 corn    rio brilhante          115.    24.8
 6 corn    maracaju               113.    22.1
 7 corn    antônio joão           110.    10.7
 8 corn    dourados               108.    38.8
 9 corn    naviraí                108.    20.7
10 corn    deodápolis              99.9   23.4
# ... with 22 more rows
g2 <- ggplot(crop, aes(x = productivity, fill = species)) +
  geom_density(alpha = 0.5)
g2

summary(crop.sum$avgProd)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  34.07   55.94   61.65   75.30  107.88  128.45 
g3 <- ggplot(data = crop.sum, aes(x = site, y = avgProd, color = species)) +
  geom_point(size = 2)
g3 + theme(axis.text.x = element_text(angle = 90, hjust = 1))

Dos gráficos, nota-se que há variação significativa entre os lugares, o ano de plantio e o tipo de cultura. De alguma maneira, todas essas variáveis deverão, então, constar do modelo final.

A seguir, unem-se os dados climáticos com os de produtividade, a fim de formar o dataset que será usado na etapa de modelagem:

clima <- read_csv("clima.csv")
clima$harvest <- factor(clima$harvest)
crop.df <- merge(crop, clima, by = 'harvest')

crop.df$material <- str_replace_all(crop.df$material, "[^[:alnum:]]", "")
seeds <- substr(crop.df$material, 1, 2)
tbl.seeds <- table(seeds)
sum(tbl.seeds/sum(tbl.seeds)*100) # check sum=100
[1] 100
barplot(tbl.seeds/sum(tbl.seeds)*100) # types of materials

crop.df$site <- str_replace_all(crop.df$site, "[^[:alnum:]]", "")
crop.df$site <- str_replace_all(crop.df$site, "[[:punct:]]", "")
crop.df$material <- seeds
crop.df <- crop.df %>% select(material, site, species, contains('rain_q'), 
                              productivity)

write_csv(crop.df, 'crop_df.csv')

O arquivo final, de nome crop_df.csv, contém os dados necessários à etapa de modelagem matemática da produtividade.

DT::datatable(crop.df)

3. Modelo preditivo da produtividade das culturas de soja e milho

library(caret)
library(ranger)
crop_df <- read_csv("crop_df.csv")
crop_df$material <- substr(crop_df$material, 1, 1) # try to balance things out
barplot(table(crop_df$material)) # still not balanced, but let's give it a shot

crop_df[, 1:14] <- lapply(crop_df[, 1:14], as.factor)

## Not run (takes a lot of time):
## Random Forest
# ctrl <- trainControl(method = "cv", selectionFunction = 'oneSE')
# grid <- expand.grid(.mtry = c(6, 8, 10), .splitrule = c("extratrees"),
#                     .min.node.size = c(1, 3, 5))
# set.seed(300)
# dim(crop_df)
# m.rf <- train(productivity ~ ., data = crop_df, method = "ranger",
#               metric = "RMSE", trControl = ctrl, tuneGrid = grid)
# m.rf
# yhat.rf <- predict(m.rf, crop_df)
# plot(crop_df$productivity, yhat.rf)
# abline(lm(yhat.rf ~ crop_df$productivity), col = 'darkred', lwd = 2)
# cor(crop_df$productivity, yhat.rf)^2
# saveRDS(m.rf, 'm_rf.rds')
## End (Not run)

# In order to save the time required to train the RF model
m.rf <- readRDS('m_rf.rds')
m.rf
Random Forest 

14436 samples
   14 predictor

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 12993, 12993, 12992, 12993, 12991, 12994, ... 
Resampling results across tuning parameters:

  mtry  min.node.size  RMSE      Rsquared   MAE      
   6    1              14.71343  0.8417227  10.933627
   6    3              14.66929  0.8426394  10.907607
   6    5              14.63400  0.8432560  10.872263
   8    1              13.18050  0.8636927   9.683304
   8    3              13.18721  0.8637343   9.695281
   8    5              13.20974  0.8634428   9.710814
  10    1              12.51664  0.8729141   9.101445
  10    3              12.50263  0.8731346   9.094314
  10    5              12.54313  0.8725464   9.126912

Tuning parameter 'splitrule' was held constant at a value of extratrees
RMSE was used to select the optimal model using  the one SE rule.
The final values used for the model were mtry = 10, splitrule =
 extratrees and min.node.size = 1.
yhat.rf <- predict(m.rf, crop_df)
plot(crop_df$productivity, yhat.rf,
     main = "Random Forest",
     xlab = "Produtividade, sacas / ha",
     ylab = "Produtividade calculada, sacas / ha")
abline(lm(yhat.rf ~ crop_df$productivity), col = 'darkred', lwd = 2)
r2 <- cor(crop_df$productivity, yhat.rf)^2
text(25, 140, paste0("R2 = ", round(r2, 3)))